#!/usr/bin/env python

# Demonstrate how to use the vtkBoxWidget to control volume rendering
# within the interior of the widget.

import vtk
from vtk.util.misc import vtkGetDataRoot
VTK_DATA_ROOT = vtkGetDataRoot()

# Load a volume, use the widget to control what's volume
# rendered. Basically the idea is that the vtkBoxWidget provides a box
# which clips the volume rendering.
v16 = vtk.vtkVolume16Reader()
v16.SetDataDimensions(256, 256)
v16.GetOutput().SetOrigin(0.0, 0.0, 0.0)
v16.SetDataByteOrderToLittleEndian()
v16.SetFilePrefix("/media/DOC/Medical-Imaging-Doc/3d/3D/DICOM/dong/MRbrain/MRbrain")
v16.SetImageRange(1, 93)
v16.SetDataSpacing(1, 1, 2)
"""
extract=vtk.vtkExtractVOI()
extract.SetInputConnection(v16.GetOutputPort())
extract.SetVOI( 0, 29, 0, 29, 0, 15)
extract.SetSampleRate(2,2,2)
v16=extract
"""
srange = v16.GetOutput().GetScalarRange()
min = srange[0]
max = srange[1]
diff = max-min
slope = 4096/diff
inter = -slope*min
shift = inter/slope
print min,max,slope,inter,shift

shifter = vtk.vtkImageShiftScale()
shifter.SetShift(shift)
shifter.SetScale(slope)
shifter.SetOutputScalarTypeToUnsignedShort()
shifter.SetInputConnection(v16.GetOutputPort())
shifter.ReleaseDataFlagOff()
shifter.Update()
v16=shifter
print v16.GetOutput().GetScalarRange()
tfun = vtk.vtkPiecewiseFunction()

tfun.AddPoint(700, 0.0)
tfun.AddPoint(20000, .0)
tfun.AddPoint(40950, 0.7)
"""
tfun.AddPoint(70.0, 0)
tfun.AddPoint(599.0, 0.0)
tfun.AddPoint(600.0, 0.1)
tfun.AddPoint(1195.0, 0.1)
tfun.AddPoint(1200, .2)
tfun.AddPoint(1300, .3)
tfun.AddPoint(2000, .3)
tfun.AddPoint(4095.0, 0.0)
"""
ctfun = vtk.vtkColorTransferFunction()
c=2

ctfun.AddRGBPoint(43000.0/c, 0.5, 0.5, 0.5)

compositeFunction = vtk.vtkVolumeRayCastCompositeFunction()
#compositeFunction = vtk.vtkVolumeRayCastMIPFunction()
volumeMapper = vtk.vtkVolumeRayCastMapper()
#volumeMapper = vtk.vtkVolumeTextureMapper2D()
#volumeMapper=vtk.vtkVolumeShearWarpMapper()

volumeMapper.SetInputConnection(v16.GetOutputPort())
volumeMapper.SetVolumeRayCastFunction(compositeFunction)
#volumeMapper.FastClassificationOn()
#volumeMapper.SetFunctionTypeToComposite()

volumeProperty = vtk.vtkVolumeProperty()
volumeProperty.SetColor(ctfun)
volumeProperty.SetScalarOpacity(tfun)
volumeProperty.SetInterpolationTypeToLinear()
volumeProperty.ShadeOn()

newvol = vtk.vtkVolume()
newvol.SetMapper(volumeMapper)
newvol.SetProperty(volumeProperty)


outline = vtk.vtkOutlineFilter()
outline.SetInputConnection(v16.GetOutputPort())
outlineMapper = vtk.vtkPolyDataMapper()
outlineMapper.SetInputConnection(outline.GetOutputPort())
outlineActor = vtk.vtkActor()
outlineActor.SetMapper(outlineMapper)


# An isosurface, or contour value of 500 is known to correspond to the
# skin of the patient. Once generated, a vtkPolyDataNormals filter is
# is used to create normals for smooth surface shading during rendering.
# The triangle stripper is used to create triangle strips from the
# isosurface these render much faster on may systems.
skinExtractor = vtk.vtkContourFilter()
skinExtractor.SetInputConnection(v16.GetOutputPort())
skinExtractor.SetValue(0, 1150)
skinNormals = vtk.vtkPolyDataNormals()
skinNormals.SetInputConnection(skinExtractor.GetOutputPort())
skinNormals.SetFeatureAngle(60.0)
skinMapper = vtk.vtkPolyDataMapper()
skinMapper.SetInputConnection(skinNormals.GetOutputPort())
skinMapper.ScalarVisibilityOff()
skin = vtk.vtkActor()
skin.SetMapper(skinMapper)


# Create the RenderWindow, Renderer and both Actors
ren = vtk.vtkRenderer()
renWin = vtk.vtkRenderWindow()
renWin.AddRenderer(ren)
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)

# The SetInteractor method is how 3D widgets are associated with the
# render window interactor. Internally, SetInteractor sets up a bunch
# of callbacks using the Command/Observer mechanism (AddObserver()).
boxWidget = vtk.vtkBoxWidget()
boxWidget.SetInteractor(iren)
boxWidget.SetPlaceFactor(1.0)

# Add the actors to the renderer, set the background and size
#ren.AddActor(outlineActor)
ren.AddVolume(newvol)
#ren.AddActor(skin)
ren.SetBackground(1, 1, 1)


  	
# When interaction starts, the requested frame rate is increased.
def StartInteraction(obj, event):
    global renWin
    renWin.SetDesiredUpdateRate(10)

# When interaction ends, the requested frame rate is decreased to
# normal levels. This causes a full resolution render to occur.
def EndInteraction(obj, event):
    global renWin
    renWin.SetDesiredUpdateRate(0.001)

# The implicit function vtkPlanes is used in conjunction with the
# volume ray cast mapper to limit which portion of the volume is
# volume rendered.
planes = vtk.vtkPlanes()
def ClipVolumeRender(obj, event):
    global planes, volumeMapper
    obj.GetPlanes(planes)
    volumeMapper.SetClippingPlanes(planes)
 

# Place the interactor initially. The output of the reader is used to
# place the box widget.
boxWidget.SetInput(v16.GetOutput())
boxWidget.PlaceWidget()
boxWidget.InsideOutOn()
boxWidget.AddObserver("StartInteractionEvent", StartInteraction)
boxWidget.AddObserver("InteractionEvent", ClipVolumeRender)
boxWidget.AddObserver("EndInteractionEvent", EndInteraction)

outlineProperty = boxWidget.GetOutlineProperty()
outlineProperty.SetRepresentationToWireframe()
outlineProperty.SetAmbient(1.0)
outlineProperty.SetAmbientColor(0, 0, 1)
outlineProperty.SetLineWidth(1)

selectedOutlineProperty = boxWidget.GetSelectedOutlineProperty()
selectedOutlineProperty.SetRepresentationToWireframe()
selectedOutlineProperty.SetAmbient(1.0)
selectedOutlineProperty.SetAmbientColor(1, 0, 0)
selectedOutlineProperty.SetLineWidth(3)


iren.Initialize()
renWin.SetWindowName("OpenGL")

renWin.Render()
iren.Start()
